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ABSTRACT 


A 336 degree of freedom slave finite element possessing capability to 
analyze engine structures under severe thermo-mechanical loading is presented. 
Description of the theoretical development and demonstration of that element 
is : presented in Volume I of this report, while a description of the 
development and use of the associated computer software is presented in Volume 
II. 
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1.0 INTRODUCTION 


This final report describes the general theory for developing slave 
finite elements, and the procedures followed for formulating element 
matrices and vectors associated with a 336 degree of freedom (d.o.f.) ele- 
ment, and the demonstration of that element, as developed under NASA 
Contract No. NAS3-23279. Description of the development and use of the 
computer software that models the capability of the element are found in 
an accompanying volume [l.l] . 

The objective of the program was to develop finite elements for use 
in nonlinear analysis of engine structures that were superior to those 
already being used for that purpose in terms of accuracy and efficiency. 

It was believed that this could be accomplished through the ’'slave 11 con- 
cept, i.e., the maximum exact use of the chosen displacement shape 
functions in the calculation of strain and stress, and the integration 
over the element volume and surfaces necessary to form element matrices 
and vectors. In addition, time was to be embedded directly into the 
formulation. This was done by formulating four dimensional space-time 
elements [1.2, 1.3, 1.4] and developing solution algorithms for their 
use [1.5] . 

In the following section, the general theory upon which all slave 
finite elements would be based is presented. The subsequent sections 
then describe the application of this theory to the specific element 
developed, and the demonstration of the performance of that element. 

The report then concludes with an assessment of the study and recommenda- 


tions for the future research. 
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2.0 SLAVE FINITE ELEMENT CONCEPT: GENERAL NOTIONS 

In this section, a detailed description of the theoretical considerations 
that would be applied to all slave finite elements will be given. In 
particular, the non-linear analysis algorithm will be described. 

2*1 Space-Time Discretizations 

The slave finite element method attempts to incorporate time as a 
dimension on a par with the spatial dimensions. It is thus reasonable to refer 
to four dimensional domains of analysis referred to as structural histories. 
Generally, one is interested in the response of the entire structure from 
some initial time, t*=0, to some final time, t=t^. Geometrically, this type 
of structural history can be thought of as a four-dimensional prism with the 
longitudinal direction being the time dimension and the three-dimensional 
cross-section being the geometric shape of the structure. 

The structure is now discretized into four dimensional finite elements in 
much the same way that conventional finite elements are formed from a 
structure. Several discretizing schemes may be used. One scheme, referred to 
as a prismatic discretization, slices the structural history at various time 
locations; then, each cross-section is divided the same way into conventional 
finite elements. An element, then is a bundle of four-dimensional longitu- 
dinal fibers extending in length 11 from t=t^ to t-t^ 9 and whose cross-section 
is the three-dimensional spatial finite element. The physical interpretation 
of this element is clear; it is a particular region in space over a particular 
time interval, and thus strongly resembles the entire domain. The element 
derived in section 3 is a prismatic element, and it is anticipated that most 
elements to be derived in future research will be prismatic. Figure 2.1 (a) 
illustrates a typical prismatic grid for bar elements, which are the two- 




(a) Prismatic Grid 


x 


(b) Semi -prismatic Grid 



Figure 2.1. Typical Grids for Space-Time Elements 
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dimensional (one spatial, one temporal) analog of a general structure. 

Other easily interpretable discretizations can be formulated 
rather easily. Figure 2.1(b) illustrates a semi-prismatic discretization. 

Note that each element is prismatic, though the space-time grid is not. 

The advantages of such a discretization would be to have more detail in 
critical space-time regions where greater accuracy is necessary, and 
correspondingly less detail in regions where the gradients in the 
response variables are relatively small. 

Figure 2.1(c) has the nodal points arranged in a prismatic manner, 
though the triangular elements shown require some thought as to their 
physical interpretation; in particular, the nature and behavior of the 
diagonal boundary line requires subtle understanding of the theory. 

The theory presented below will assume that the elements are 
prismatic in nature, even if the discretization is not. It is believed 
that this offers facility in both physical interpretation and element 
formulation. 

2.2 Hamilton's Law of Varying Action 

The basic principle upon which slave finite element matrices will 
be formulated is Hamilton* s law of varying action [2.1]. In order to 
derive this law, one begins with the appropriate form of the variational 
statement of the principle of virtual work for the particular structural 
member or element: 

f a • i £ dV = ; *£*• 6 u dV + / T • 6 u dS (. 

v v sa 



where u, e, a, and T represent the generalized displacement, strain, 
stress and traction fields appropriate for the element. The term f* 
Includes conventional and inertial body forces, and can be expressed as 


f* - f - (2.2) 

where p is the momentum per unit volume of a material point in the body, 
measured relative to an inertial reference frame, in excess of the rigid 
body motion of the structure. The inertial loads due to the rigid body 
motion are incorporated into f. In addition, it should be noted that 
the time derivative is taken in the inertial frame. 

Substituting (2.2) into (2.1) and bringing the momentum term to 
the left-hand side of the equation yields 


/ (o*6e + 4?- • 6 u) dV * / f • 6 u dV + / T • 6 u dS 
V dt V Sc 


(2.3) 


Both sides of (2.3) are integrated with respect to time from t=t^ to 
t-t^. Noting that the velocity V, relative to the inertial reference frame 
in excess of rigid body motion, is defined as 


du 

dt 


(2.4) 


the momentum term may be integrated by parts over time to yield 

c 2 — _ _ _ C 2 t 2 

/ -r^~ *6udt = p*6u | - / 'p'Sv dt 

t dt t t 

1 Z 1 Z 1 


(2.5) 


Hamilton's law of varying action is thus stated as 


f J (o*6e-p*6v) dVdt- = f f f*6u dVdt + / / T*6u dSdt - / p*6u|dV 
t x V V t 1 So V t x 


( 2 . 6 ) 
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The last two terms on the right side of (2.6) are actually similar, in that 
they represent the work done 'by the "loads" on the three-dimensional 
hypersurface of the element. The first of these terms covers the surfaces 
whose outward normal are perpendicular to the time dimension, while the 
second of these terms cover the "prism ends", or, what will be referred 
to as the time boundary. Note that if u(t2) or u(t^) are specified, that 
one or both of the terms in this integral would be zero. 

2.3 Typical Linear Elastic Element Formulation 

While the linear elastic case for any of the slave finite elements 
must be derived as a special case of an element normally used in non-linear 
analysis, it is useful to study a straight-forward derivation in the 
elastic range in order to grasp some of the basic concepts of the slave 
finite element method. 

Let u be expressible as 

u = [n] ( u) (2.7) 

where [N] are a set of shape functions and {u} are degrees of freedom (d.o.f.) 
of the element, generally nodal displacements. Now, unlike conventional 
finite element transient analysis, where an individual element of {u} 
represents a function of time indicating the response at a particular point 
in the structure, an element of {u} represents the response of the structure 
at a particular point in space and time; thus, {u} is really a set of 
constant coefficients to be determined only by algebraic equations. 

The strain and velocity fields are defined as 

t = [N']{u} (2.8a) 

V = [N ]{u} 


(2.8b) 
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where [N - *] and [N] are the shape functions [N] acted on by assumed linear, 
first order, differential operators. Since these quantities will represent 
the highest order derivatives of u to be used in the left-hand side of 
(2.6), the shape functions [N] must satisfy interelement continuity. 

The linear constitutive laws are 


0 - [E] e (2.9a) 

P= [p] V (2.9b) 

where [e] is referred to as the modulus matrix and [p] as the mass density 
matrix. Using (2.8) and (2.9) in the left-hand side of (2.6), and taking 
the variations on {u} yields 

C 2 _*_*_*_» T 

/ / (0*5e - p»6v) dVdt = {6u} [K] (u) (2.10) 

Z 1 V 


where the dynamic stiffness matrix [K] is expressible as 


[K] = / 
t. 


/ ([m T [E][in - [N] T [p][N]) dVdt 
V 


( 2 . 11 ) 


Note that for [E] and [p] symmetric that [K] is symmetric. Also, 
assuming [E] and [p] are positive definite that [K] is the difference 
between two positive definite matrices. Calling the first term the 
stiffness contribution and the second the inertial contribution, it is 
found that, generally, the stiffness contribution increases with 
for a fixed volume, while the inertial contribution decreases. This 
divides the types of problems qualitatively into three regions: 
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i) T "large" - In this case, the stiffness terms totally 
overwhelm the inertial terms, which may be effectively 
ignored. This case, then, is appropriate for quasi- 
static analysis. Many problems involving engine structures 
would fall into this category. 

ii) T "small" - In this case, the inertial terms dominate the 
analysis. Certainly, a series of rapid pulses, with a 
period significantly less than the fundamental period of 
the structure, -may be analyzed using a dynamic stiffness 
composed only of inertial terms. 

iii) T Moderate 11 - With the stiffness and inertial terms 

comparable, effective vibrational analysis can be performed. 
Depending on the sophistication of the shape functions, 
particularly in time, values of T on the order of .10 to 
.25 periods are not unreasonable. (To determine the 
length of a period, see Section 2.8 below.) 

The choice of whether to use i) , ii) or iii) is dependent upon the time 
variation of the loading relative to the period of vibration of the 
fundamental mode of the structure. 
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The right-hand side of (2.6) can be collected into a force vector 
of which some terms are known and others are not. For instance, the body 
force term, which includes the rigid body inertial loading, may be 
transformed into a known vector by 



/ f“*6u dVdt 
V 


{6u} T (F b } 


where 


{ V 


/ 

t. 


f [N] T t dVdt 
V 


Generally, for elements internal to the space-time domain, the other 
two terms on the right-hand side of (2.6) represent tractions and/or 
impulses per unit volume that neighboring elements exert on the element 
under study. Their vector equivalent is defined as the (unknown) {F}. 
Combining all such {F} vectors into {F*} yields the element equations 
[K] (u) = {F*} 

Note that if all the\u^ have dimensions of length that {F*} has units of 
force times time. The physical interpretation of {F*} is that they are 
the equivalent point impulse-momentum difference to the actual distributed 
loads placed on the element. 

The equation (2.14) resembles the equation of statics in conventional 
finite elements. It is well known that the stiffness matrix for a con- 
ventional finite element is singular. The rank of the matrix is reduced 
by the number of equilibrium equations available to the element. It 
will now be shown that similar conditions exist for the system (2.14). 


(2.12 


(2.13 


(2.14 
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The differential equations of motion are 
°ij t j + f i = 

Integrate this expression over the space-time domain. Noting that the 
tractions are related to the stresses on a surface whose outward normal 
is n by 

T i ■ °u "j 

and employing the divergence therem, one obtains for the left side of 
the equation 


I . = / 4 T, dSdt + / / f dVdt 
1 4 s a 1 t x V 

where it is seen that 1^ represents the total impulse over the element. 
The right-hand side is shown to be 


AP = / p. | dV 

V C 1 

or, the change in linear momentum. Thus, the linear impulse-momentum 
equations must be obeyed. Similarly, by taking moments, 

tf2 l *iai x i f i dVdt+ / 2 t £ kii x i T i dsdt = f v e kii x i p i U' dv 


c i v 


c i s 


or, that the angular impulse -moment urn equations must be satisfied. 

These two sets reduce the rank of the matrix by 6. 

If the inertial effects are not important, then the impulse-momentum 
equations reduce to 


f f f dVdt + / i T dSdt = 0 
ti V 1 t x s 


1 1 e kli X 1 f i dvdt + ; 6 c kli X 1 T i dsdt = 0 

C 1 v t l s 


(2.15) 


(2.16) 


(2.17) 


(2.18! 


(2.19. 


(2.2a 


( 2.20 
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Note that in (2.20a) and (2.20b) a sum of spatial integrals can be brought 
under one time integral. In order to be true for arbitrary time histories, 
the sum of the spatial integrals must be zero for all times. Physically, 
this means that the body is alvays in static equilibrium, which, of course, 
is exactly what is desired in the quasi-static problem. Furthermore, 
this may reduce the rank of [K] by 6 for each time station. 


2.4 Adjustment for Non-Interelement Continuous Shape Functions 

Many times, a set of shape functions [N] satisfying inter-element 

continuity may not be easily obtainable. One method which can help 

maintain continuity would introduce into the formulation a displacement 

field Ug defined on S u which would satisfy the continuity requirements 

once specification of the (u) was made. The displacement field u would 

then be constrained to Ug on (in the average sense), with the tractions 

on serving as weighting functions; thus an additional term 

t 2 _ __ 

6 (/ / T* (u-Ug)dSdt) 

t l S u 

is summarized from the left-hand side of (2.6). The boundary displacement 
Ug is expressed as 

"b = [n b ]{ u} 

and the boundary functions T are expressed in terms of {u} through the 
strain-displacement and constitutive laws evaluated on S^. The result is 
that (2.21) is quadratic in {u}, and, after variations are taken, a 
symmetric matrix, composed of the sum of 

[x] = / 2 / [N'] T [E]([Ng]-[N])dSdt 


t. s 

1 u 


and its transpose is added to [K] to form the stiffness matrix for the 


( 2 . 21 ) 


( 2 . 22 ) 


(2.23) 


element. 
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2.5 Global Assembly Considerations 

For most problems, assembly of equations (2.14) into a global system 
Is fairly straight-forward. The same principle used to assemble a global 
stiffness matrix in ordinary finite elements can be used here, i.e., 
identification of a local d.o.f. with a global d.o.f. The right-hand 
side of the equation is developed by calculating the contributions due 
to body forces, tractions, line and point loads, impulses per unit volumes, 
and point impulses into a vector {P}; thus, 


EWiH“> ■ {pJ 

Boundary conditions (both spatial and temporal) reduce the system (2.24) 
to the one that models the structural history, and then the system is 
solved. 

The spatial boundary conditions generally pose no problems. If 
the spatial boundary conditions are posed properly from a solid mechanics 
viewpoint, that is, + S u = S for each conjugate pair that could be 


specified, then elimination of corresponding rows and columns of [K 
is performed in the same manner as ordinary finite element analysis. 


global 


] 


The temporal boundary conditions consistent with this formulation 
would be the specification of the displacement or the momentum at t=0 and 
t=tp. In order to uniquely determine the solution, the displacement must 
be specified at least at one of these times. Unfortunately, the classical 
transient analysis problem is usually specified in terms of initial 
conditions , which is the equivalent of specifying the displacement and 
the momentum at t=0, and not specifying either at t=tp. The approach to 
properly incorporating initial conditions depends on the temporal shape 


(2.24) 


functions used. 
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2.5 >1 Lagrangian Elements 

By Lagrangian elements it is meant that only displacements and not 
any of their (time) derivatives are specified as d.o.f. Specifying initial 
conditions would then be equivalent to specifying a boundary condition 
and the corresponding reaction force, without specifying an active 
force at the other end of the one structure. 

The time boundary terms are volume integrals of the momentum and 
variations of displacement at t=0 and t=tp. Assume the initial conditions 
are in form 

u (o ) - u q = 0 (2.25a) 

p(o) - P q = 0 (2.25b) 

Condition (2.25b) is incorporated into the time boundary term with a 
language multiplier field X as 


/ (p(t )*6u(t ) - p(o) *6u(o) + 5[A«(p(o) - P )]) dV 
V - r r o 


Variations are taken with respect to X and p(o) to yield 


(2.26) 


/ (p(t F )*6u(t ) - p(o) *6u(o) + (p(o)-p ) -6X +X*6p(o) ) dV (2.27) 

v r r ° 

If X is identified as u(o)-u o , then (2.27) becomes 


/ (p(t F ) •6u(t r ) - p 6u(o) + (u(o)-u )6p(o)) dV 
v F F o o 


(2.28) 


The implications of (2.28) are that u(o) is treated as an independent variable 
satisfying (2.25a); that this constraint is satisfied in a straight forward 
way so that it may be easily substituted in (2.24), thus eliminating the columns 
of corresponding to u(o) ; on the other hand, since it is assumed 
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that Su(o)i*0, the rows corresponding to P(o) are not eliminated, but are 
maintained with p Q incorporated properly into the right-hand side of (2.24). 

-V. 

Replacing u(o) as unknown in the problem is P(tp). 

There are a few pitfalls in this approach. First, the number of d.o.f. 
describing u(o) and P(t_) must be the same. Secondly, the symmetry of 
[Kgiobai^ * s not maintained. Thirdly, since the equations at t»0 are 
maintained in the analysis, a prismatic discretization scheme results in 
equations equivalent to explicit time integration schemes. Some non-prismatic 
discretizations may result in no solution being obtainable. Some of these 
problems may be solved by either additional constraints or solving in the "least 
squares" sense, but many would find this distasteful. 

2.5.2 Hermitian Elements 

By Hermitian elements, it is meant that both displacements and velocities 
are specified as d.o.f. These elements are at least cubic in time, and thus 
are quite accurate. Through the constitutive law (2.9b) » it is seen that sped- 
fying p is equivalent to specifying v; thus, certainly, initial conditions would 
be easily specified; however, several changes must be made. Firstly, the 
time boundary terms no longer contribute to the right-hand side of the equation. 
Instead, ^ K g^ 0 bal^ is ac *j ustec * by adding to the rows corresponding to the 
displacements at the t=tp and the columns corresponding to the velocities at 
t=tp the conventional finite element consistent mass matrix; and subtracting 
this matrix from the rows corresponding to the displacements at t=0 and the 
columns corresponding to the velocities at t*0. 
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The rows and columns corresponding to u(o) are eliminated by 
specifying the initial displacements. This seems reasonable and proper. 

It also eliminates one of the off-diagonal adjustments noted in the 
paragraph above. To incorporate the initial velocity condition, the 
arguments of [2.1] will be followed. 

The effect of the remaining off-diagonal adjustments are to cause a 
degeneracy in C K gi 0 bal^ a ^ ter t * ie initial displacement condition is 
imposed. How well or poorly this is accomplished will affect the behavior 
of the solution. Ideally, elimination of the rows and columns corres- 
ponding to V(o) would be the best choice from a physical viewpoint, 
because of the specification of V(o), and from the computational viewpoint, 
since [K would remain nearly symmetric. Unfortunately, some simple 

free vibration problems exhibit anamolous damping with refinement in the 
grid spacing. This is considered due to an inability of the additional off- 
diagonal term to render exactly singular. The situation is 

remedied by asserting that the rows corresponding to the final displace- 
ments should be linearly dependent on all the remaining rows of the matrix. 
After asserting that 6u(tp)^0, and substituting a linear combination of 
the remaining rows for the rows corresponding to the 6u(tp) equations, it 

is found that the non-trivial portion of the 6u(t-) equations are those 

r 

equations associated with 6V(o) . Thus, the rows associated with the 


6u(tp) equations and the columns associated with V(o) are eliminated, 

rendering [K in its final form. Note that the symmetrical nature 

of EKg]_ 0 b a ]_3 is destroyed. Note further that this procedure works if the 

number of variables in V(o) and u(t ) is the same, implying that prismatic 

r 

discretizations would be preferred. On the other hand, the remaining 
equations do not form an explicit solution scheme pattern^ lending stability 


to the process. 
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The other pitfalls of Hermitian elements include the fact that 
velocities are indeed discontinuous in the presence of impulsive loading. 
This difficulty is handled by introducing 2 sets of velocities, one 
slightly before and one slightly after impulse, and relating them by a 
work-equivalent constraint involving the impulse-momentum equations. In 
addition, while the point-loads conjugate to the displacement d.o.f. 
can be thought of as pure impulses, the "time moments" conjugate to 
velocity have no real physical meaning. 

2.6 Non-linear Analysis 

The discussion presented thus far acquaints the reader with basic 
notions about the slave finite element concept, and discusses the techniques 
of linear anlaysis. It should be mentioned that inversion of the final 
form of (2.24), after boundary and initial conditions are accounted for, 
results in the solution of the entire transient problem . 

2.6.1 General Approach 

Like most non-linear algorithms, the one presented here is based on 
an iterative procedure involving quasi-linearization. The difference in 
philosophy between the method presented here and the conventional step-by- 
step methods, such as the tangent modulus or residual force methods, is 
demonstrated graphically in Figure 2.2. In both Figures 2.2(a) and 2.2(b), 
the heavily drawn curve represents the "exact" time history of quantity A, 
which may be thought of as the displacement of a certain point, or a 
stress at a point, etc. In Figure 2.2(a), the step-by-step method has 
calculated the path OA, representing the history up to that point. At 
A, the procedure generates successive approximations AB^ until convergencw 
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to the path AB is achieved. The procedure is now repeated at B. For 
the algorithm presented here. Figure 2.2(b) indicates that successive 
iterations generate an entire load history. This load history allows for 
determining loading and unloading paths for the entire time interval 
of interest, eliminating the guesswork or sub-interval changes usually 
associated with the step-by-step methods. Iterations converge until the 
curve OC is obtained. 

2.6.2 Element Matrices 

The derivation of element matrices will be based upon proper 
substitution into (2.6). The displacement field is expressed as in 
(2.7) and is assumed admissable. Finally, the time interval is trans- 
lated so that 0<t<T. 

The strain-displacement law may be written formally as 
e = f (u) 

where f is a function of u and its spatial derivatives, and, in general, 
is non-linear. Taking variations of (2.29) yields 

5e - f"(u)6u 

where the prime has a general meaning relating to derivatives with respect 
to u and its partial derivatives. In a similar manner, the velocity- 
displacement relations may be expressed as 

V = g(u) 

where g, like f is a (non-linear) function of u and its spatial and 
temporal first derivatives. Taking the variation of (2.31) yeilds 

-Jh. _Jk. 

6v = g'(u)6u 


(2.29) 


(2.30) 


(2.31) 


(2.32) 


It is also interesting to note the time derivatives of (2.29) and (2.31): 
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o . . o 

— » s X x ^ 

e f (u)u 

V g (u)u 

It will be assumed that the constitutive laws will be linear between 
the time derivatives of the appropriate kinematical and dynamical 
quantities. Some linearization must be available in order for the matrix 
techniques of finite element analysis to be applicable. As a result, 

the laws relating to e and o, and v and p are formulated as 

5. _ 2.+ 3. _ _ _ _ 

0 [E(0, C, 0, X, t , . . .)3 ^ T (Q t E , 0, X, t , . « . ) 

Si = o. + £ 

P [p(0, E > 9» t,.. )]v ^ ( 0 * E , 6 , x, t , . • . ) 


(2.33a) 

(2.33b) 


(2.34a) 

(2.34b) 


where 0 is the temperature field. Equations (2,34) are integrated from 
0 to t; thus. 


t £ - 

a = f o dt + o 

0 

t _o ^ ^ 

p = / pdt + p 

A U 


(2.35a) 

(2.35b) 


where a and p are the values of the stress and momentum at local time 
o r o 

t=0. These may be approximated by 

O o “ [N o (x)](a o } (2.36a) 

p o = [N p (x)]{p Q } (2.36b) 

where the number of parameters in {a Q } and {p Q } is aribtrary. Ideally, [N ] 
is inter-element continuous and both [N^] and [N^]are as sophisticated as 
the stress and momentum fields derived from [N] in the linear theory, 
though neither requirement is really mandatory. 
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In the iterative procedure to be used, "current" values of 

displacement, stress, strain* etc., for the entire load history are 

in hand. These values are used to evaluate i' , g' , [E], t, [p] and 
o o 

it. The quantities 6u and u are calculated from (2.7) where the set {u} 
are unknown, representing the "updated" solution. Similarly, {a Q } and 
{p Q } are unknowns. 

Equations (2.30), (2.32) and (2.33) may be represented as 
fie - [f'][N]{fiu} 

e = [f'][N](u) 

<5v = [g'][N]{Su} 
v = [g'][N]{u} 

where [f'] and [g'] are 6x3 and 3x3 operator matrices, respectively. 
Stress and momentum matrices are defined by 

[S] = /[E][f'][N]dt' 

0 

[M] = / [p][g-][N]dt- 

0 

and stress and momentum vectors are defined as 
t c> 

T - / Tdt' 

0 

-* = c 
n f v dt' 

0 

Using (2.36), (2.38) and (2.39) in (2.35) yields 


(2.37a) 

(2.37b) 

(2.37c) 

(2.37d) 

(2.38a) 

(2.38b) 

(2.39a) 

(2.39b) 
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a = [S]{u} + T + [N a ](a o } (2.40a) 

p = [M]{u> + tT +.[N p ]{p o ) (2.40b) 

Equations (2.4) may be evaluated at t=T. (Quantities evaluated 
at this time are given a T subscript.) Assuming that the shape functions 
[N 0 ] and [N ] can also approximate 0^ and p^, then equations (2.40) take 
on the form 

[S T ](u} + t“ t + [N 0 ](a o - a T ) = {0} (2.41a) 

[Hj.Hu} + }T t + [N p ](p o - p T ) = {0} (2.41b) 

at t=T. Equations (2.41) are used as subsidiary conditions to the problem. 


They are added to (2.6) with the use of languangian multiplier fields. In 
particular, the choices for the stress and momentum fields, respectively, 
are given as 


X 

P 


■ cw 

- rv { V 


(2.42a) 

(2.42b) 


These fields are used with (2.41) and integrated over the volume of the 
element. These conditions, as well as (2.37a), (2.37c) and (2.40) are 
used in (2.6) to give 

(6u} T ([K e ]{u} - {F*} - {F*} - {ip + [Ae]{a o } - [A*]{p o }) 

+ 6(tt a } T ([ B*]{u> + [C^](a o - o T >-{qp) 

+ {X p ) T ([B p ] (u) + [C p ] (p Q - P T )-{q^})) = 0 


(2.43) 
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where 


11 

1 — 1 
0) 

1 1 

/ / ([N] T [f'] T [s] - [N] T Cg'] T [M])dVdt 

0 V 

(2.44a) 

u 

1 — 1 
Q > O 
< 
i i 

/ / [N] T [f ‘*] T [N ]dVdt 
0 V 

(2.44b) 

i — 1 

1 1 

II 

T 

/ f [N] T [g'] T [N ]dVdt 
0 V p 

(2.44c) 

ttf- 

: [n ] T [s ]dv 
v 

(2.44d) 

1 — 1 
w 

*T3 fD 

\ I 

tl 

/ [Np] T [M T ]dV 

(2.44e) 


f [N ] T [N]dV 
v a a 

(2.44f) 

[C P- 

f [N ] T [N ]dV 
v p p 

(2.44g) 

< F t» ' 

T 

-I I [N] T [r] T TdVdt 
0 V 

(2.44h) 


T 

S / [N] T [g“’j r TTdVdt 

0 v 

(2.44i) 

<S J ■ 

-/ [N ] t t T dV 
V 

(2.44j; 

«c> - 

-f [N ] T tt dV 
v P T 

(2.44k; 


The set of equations generated by (2.43) when variations are taken on 

{u},{a }, {p }, (X „} and (X } are 
o o o p 

[K e ](o) + [Ap{0 0 } - [A*]{p o ) + tB']\> + [Bp T U p > 

= {F*} + {F®} + {F®} 

T TT 


(2 .45a 
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rc*]u 0 ) - (o) 


(2.45b) 

[cj]ftp) - (0) 


(2.45c) 

[#{“> + CcJ]{o o 

- 0 T ) - <q*) 

(2.45d) 

[Bp{u) + [C*]{p o 

- v ■ <<> 

(2.45e) 



the multipliers identically zero. They may then be omitted from (2.45a). 

If [N] is not interelement continuous, equations (2.45a), (2.45d) 
and (2.45e) may be used (without the multipliers), with the following 
adjustments : 


EKJ = (expression (2.44a)) 4- [X g ] + [X g ] -[B®] [C Q ] [a“] 

T 

{F e } = (expression (2.44h)) + f f [N-N n ] T [r] T dSdt 

0 s u B 

[A*] = (expression (2.44h)) + [A^] 

where 

T 

[XI - / / [N -N] T [r][s] dSdt 

e o S u B 

[a"] = / / [N B -N] T [r][N n ] dSdt 

O o s B 0 


(2.46a) 

(2.46b) 

(2.46c) 

(2.47a) 


(2.47b) 



25 


2.6.3 Solution Algorithm 

Equations (2.45 a, d, e) are assembled into the systems 
[K]{u} + [A a ](a} - [A p ]{p} - {p} 

[B a ]{u} + [C - (q T ) 

[B p ]{u} + [C p ]{p} « {q^} 

where the assembly enforces the conditions of continuity across a time 
boundary for stress and continuity of momentum across a time boundary to 
the extent that no impulses per unit volume are applied, and, if such 
impulses areq>plied, an appropriate discontinuity is maintained. The 
results should leave [C^J and [C p ] square and invertible. Thus, equations 
(2.48b, c) are solved for (c) and {p} and used in (2.48a) to derive the 
global stiffness equations 
[K]{u) = {P} 

where 

[K] = [K] - [A a ][C a ] _1 [B o ] + [A p ][C p ] _1 [B p ] 

(p) = (p) - [A a ][c a ]- 1 (q T } + [A p ][ Cp r\y 

Equations (2.49) are solved for {u} and then back-substituted into all the 
pertinent equations to calculate the important quantities to be used in 
the next iteration. In particular, {a} and (p } are found from (2.48b,c). 

On an element level, the displacement u, the strain e, the velocity v, 

Ok 

the stress O and the momentum p are found from (2.7), (2.29), (2.31), 

(2.40a) and (2.40b), respectively. These new values are now used in (2.34) 

o o 

to determine new values for [E], T, [p] and ir, and in (2.37) to determine 


(2.48a) 

(2.48b) 

(2.48c) 


(2.49) 

(2.50a) 

(2.50b) 


values for operators [f'’] and Cg'*]. The process is repeated until some 
measure of convergence on the structural history is met. 
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2.7 Problems of Free Vibration and Buckling 

The problem of free vibration has meaning only in the linear range; thus 
the techniques of §2.3 are applicable. 

In conventional finite element analysis, the problem of free vibration 
reduces, to a set of linear, homogeneous, ordinary differential equations 
after the boundary conditions have been imposed. A periodic solution with a para- 
meter representing the frequency is assumed. This, of course, will yield linear, 
homogeneous, algebraic equations for the magnitudes of the d.o.f. (thus, the 
mode shape), which has a non-trivial solution only for certain values of the 
frequency parameter. If there are N d.o.f. unknowns, then there in general will 
be N such values, serving as approximations for the first N frequencies of the 
structure. Associated with each of the N frequencies is a particular vibration 
mode of the structure. 

In slave finite element analysis, the global equations (2.24) are already in 
algebraic form. It is thus necessary to impose the condition of periodicity in 
same way consistent with the formulation. This is done by stating that u(o)=u(tp)=0; 
that is, the clock begins ticking for a structure in free vibration when the 
deflections are zero; and that at some later unknown time t^,, the structure has 

r 

returned to that undeformed state; thus, the time domain’s length is the basic 

parameter of the study. It represents an integral number of 1/2 periods. The 

conditions stated above are not initial conditions, but the type of boundary 

condition appropriate for use in Hamilton’s principle. (Note that non-prismatic 

discretizations could be feasible in this case, though not necessarily desirable.) 

This renders [K square and symmetric and the right-hand side zero. Values 

for the t that render [K . , - ] singular are now found. If the grid is strictly 
r global 

prismatic, then for N spatial d.o.f. and K time station* there are NK values of tp. 
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It is found that this represents K individual approximations for the N different 
modes and frequencies. This implies that many time stations are probably not 
necessary or desirable to find approximate vibrational characteristics. 

Like free vibration problems, buckling problems reduce to a set of linear 
homogeneous equations, this time for the buckling mode. In addition, there is 
the solution for the so-called fundamental mode, that is, the form of the displace- 
ment field that occurs when no buckling is present. This mode is particularly 
important in the inelastic range, as it will determine the moduli for the con- 
stitutive laws in the buckling equations. 

The slave finite element method works well with so-called dynamic buckling 
problems [2.2]. Here, the load P(t) is expressed as 

P(t) = P o f(t) (2.51) 

where f(t) has a maximum value of 1. Note that this is different than static 
buckling, where f(t) would be monotonically increasing. Furthermore, for the right 
values of P q , the dynamic buckling phenomenon takes place immediately at t=0. 

The buckling load P q is the maximum load that allows the buckling mode to be a 
free vibration mode. Generally, it is found by setting up the buckling equations, 
with inertia included, over the interval 0<t<t P , with u(o)=u(t r )=0, t_ itself un- 
known. Since the resulting equations are linear and homogeneous, the determinant 
of the dynamic stiffness matrix is zero. This yields an expression for P q in 
terms of tp, the maximum of which is determined. In most instances, this maximum 

will occur in the limit as t + 00 . 

F 

Static buckling problems could be solved in the following manner. Generally, 
the dynamic buckling load for f(t)=l is the same as the static buckling load, 
at least in the elastic range. It is probably a reasonable assumption that the 
method extends to the plastic range, if not the result. By this, it is meant that 
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the static buckling problem is set up as a quasi-elastic dynamic buckling problem 
with f(t)=l, with the moduli determined from the fundamental solution at the 
assumed bifurcation point. On the other hand, the true dynamic buckling problem 
history runs simultaneously with the fundamental solution history, and is 
affected by it as a function of time. 

2.8 Slave Finite Element Method 

The methods described above will be used to formulate and develop slave 
finite elements. In addition to the general theoretical methods, the following 
restrictions, which comprise the basis of the "slave" approach, are placed on 
the study: 

Firstly, all quantities in the integrands of the derived matrices must be 
expressible in terms of interpolation polynomials. Generally, if the poly- 
nomials can be directly computed from manipulation on shape functions [N] (or, 

[No] or [N ]) then these polynomials will be used directly. If direct manipulation 
does not yield polynomial expressions, then values of the new quantity at the 
nodal points are taken and a new interpolation polynomial is given. 

All integrations must be exact and in "closed form." For irregular boundaries 
(non-rectangular) , this means that isoparametric representation techniques cannot 
be used; thus, the corrective matrix technique of §2.4 must be employed. 

It would also appear that prismatic elements offer the most straight- 
forward approaches to time embedment. All elements would be restricted to this 


type. 
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3.0 THE 336-D.O.F. , SLAVE FINITE ELEMENT 

In this section, a slave finite element appropriate for use 
in transient analysis for thick and moderately thin plates is 
developed in detail. The element is capable of incorporating 
plasticity and creep, as well as elasticity, and also allows for 
variation of the parameters with temperature. The structure 
may be loaded by body forces (including inertial effects), traction 
loads, point loads, body impulses and point impulses, as well 
as loads induced by thermal effects. The element will be shown 
to be compliant with §2.8 of this report. 

3.1 Element Geometry 

The element geometry is depicted in Figure 3.1. Note that 
the coordinates shown are "local" coordinates, though the time 
axis is parallel to the global time axis. The x-y plane is 
defined as the midsurface of the plate. The thickness of the 
plate is a function h(x,y), which will be limited to those 
functions having linear variation along an element edge. The 
upper and lower surfaces of the plate have projections equal to 
the midsurface area of the plate, i.e., each edge plane contains 
lines parallel to the z-axis. 

Figures 3.2 show the various boundary "surfaces" of the 
element. Figure 3.2a depicts the "surface boundary". This may 
be the upper or lower surface. It is assumed that elements of this 
type will not be connected to each other through the upper or lower 
surfaces; as a result, it serves as a traction boundary. Furthermore 
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FIGURE 3.1* 


The 336-DOF Slave Element 
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interelement continuity conditions may be relaxed somewhat on these 
surfaces. Figure 3.2b depicts the "edge boundary". Note the use 
of the edge coordinate £. Generally, the elements are connected 
spatially through the edges; however, non-rectangular geometry 
will generally lead to non-interelement continuous displacement fields. 
It is here that the corrective matrix technique of §2.4 will be 
employed. At times, however, these edges will lie on the outer 
surface of the structure and also serve as traction boundaries. 

Finally, Figure 3.2c depicts the "time boundary" which is nothing 
more than the 3-D spatial representation of the element frozen in 
time. Since the prismatic discretization leaves time perpendicular 
to space, it is likely that the shape functions will be inter-element 
continuous across a time boundary. 

The element has 56 nodal points, 28 at t=0 and 28 at t=tp. Of 
the 28 at either time station, 12 lie on the upper surface, 12 on 
the lower surface and 4 at the corners of the midsurface. This 
places the element in a four-dimensional serendipity family, as 
nodes are only on element edges. Each node has 6 degrees of 
freedom: the 3 linear displacements, u, v and w in the x, y and 

z directions, respectively, and their first derivatives with respect 
to time, denoted u, v and w. Thus, the element contains a total of 
336 dof. 


3.2 Choice of Shape Functions 

The displacement fields u, v and w are considered independent 
unknowns. As a result, each field should be expressible in terms of 
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112 independent polynomial terms in x, y, z and t. Let the collection 
of these polynomials be denoted as [p]; then the displacement field 
is expressed as 

u - [p]{a> (3.1) 

By evaluating [p] (or its time derivative, [p]) at the 56 nodal points, 
a set of relationships relating the constants {a} associated with 
the linear combination (3.1) and the set of physical dof (u) can 
be written as 

{u} = [B]{a) (3.2) 

Since the number of functions in [p] has been chosen to be equal 
to the number of dof, [B] is square, and also invertible. Solving 
for {a} in (3.2) and backsubstituting into (3.1), and noting (2.7), 
the shape functions [N] are thus 

[N] = [p][B] -1 (3.3) 

In theory, [B] is 336x336; however, since each of the displacement 
fields will be modeled using the same shape functions, [B] can 
immediately be reduced to 112x112; furthermore, because of the 
prismatic nature of the element , and the use of velocity dof, each 
shape function may be expressed as the product of one of 28 
spatial shape functions times one of 4 temporal shape functions. 

The temporal shape functions are the well known cubic Hermitian 


shape functions 

N - l-3(t/T) 2 + 2(t/T) 3 
C 1 

N - T{(t/T)-2(t/T) 2 +(t/T) 3 } 
C 2 

N = 3(t/T) 2 -2(t/T) 3 
C 3 

N = T(-(t/T) 2 +(t/T) 3 J 
C 4 


(3.4) 




The 28 spatial shape functions chosen are those that are present 
when the spatial configuration reduces to a rectangular parallelipiped. 
These shape functions are 
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Thus, only a 28x28 matrix needs to be inverted. The division of the 
shape functions into those formed by combinations in (3.5) and 
those of (3.4) guarantee that, at t=0 and t=T, the displacement field 
at any point in the body is a function only of the displacement 
dof at that time station. Similarly, the velocity field at 
any point in the body is a function only of the velocity dof at 
that time station. As a result, the kind of continuity desired 
across a time boundary is present. 

If the thickness is not constant, the displacement field on 
the upper and lower surfaces will depend on all the dof. Generally, 
this means that the displacement field will not be intereleraent 
continuous; however, since elements will not be connected to each 
other on these surfaces, no further adjustment is necessary. 
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If the edges are not parallel to the coordinate axes (x and y) , 
then the displacement field will not be interelement continuous across 
that boundary. Since elements will connect across 
this boundary, the methods of §2.4 must be used. First, note that 
[N] must be evaluated on these edges. Noting Figure 3.2b, it is 
seen that for an edge connecting corner i to corner j , that 

x s. x^ + £cos9 (3.6) 

y = y i + £sin0 


where 


e 


tan 


-1 




x.-x 

J 


i 

i 


Expressions (3.6) are substituted into [N] (or, alternatively, [p]), 
to express [N] on a given edge as functions of £, z and t. On the 
other hand, [N„] is formulated directly in terms of these variables 
by treating Figure 3.2b as if it were its own element. For the 20 
nodes showing, there are 40 dof to describe each displacement 
field. The time dependence is chosen as the same as for [N], i.e., 
(3.4); thus, 10 spatial shape functions are found by proper linear 
combinations of the terms 


i z z 2 z 3 z 

Zz Z 2 z Z 3 z z 2 Zz 2 


(3.7) 


(3.8) 


It should be noticed that on the lines £=0 and Z~Z^ this displacement 
field is interelement continuous, and thus no further adjustments 
to account for interelement continuity requirements are necessary. 
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The choices for [N^] and [Np] appear to be arbitrary. For 
convenience, they will bd derived by assuming each of the 6 stress 
fields and 3 momentum fields can be expressed using combinations of 
the 28 spatial polynomials in (3.5). This yields greater accuracy 
for the stresses than the elastic case yields and the same accuracy 
for the momentum that the elastic case yields. 


3.3 Non-Linear Strain Displacement Law 

The non-linear strain displacement laws are those determined 


by using the Green's strain tensor; i.e., 

£ ij = 2 (u i,j + U j,i + “k.j "k.i 5 
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(3.10) 


Taking variations with respect to the displacements yields 
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These equations are the particular form for equations (2.30). When 
expressed alternatively as in (2.37a), it is seen that the [f"] matrix is 
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8 , n , 9v . 

3J 3y + a 3y } 
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3x 3y 


3w 3 
3y 3x 


Sii 3 + 3ii _3_ 
3y 3z 3z 3y 




3v 3 
3z 3y 


( 1 + lr ) 


_3_ 

3y 


|H 3 + (1+ |H. ) 

3z 3x 3x 


_3_ 

3z 


3v _3 , 3_v _3_ ... 3w. 3 . 3w 3 

3z 3x 3x 3z 1 3z ; 3x 3x 3z 


(3.12) 


Note that [f'*] can be divided into strictly linear and non-linear parts. 
When both the displacements and strains are very small, the non-linear 
portion may be deleted. 
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In using (3.12) in the non-linear algorithm, it should be noted 
that the terms involving. u, v and w are found from the just completed 
iteration that solved for {u}. Using this expression back in (2.7) 
gives specific forms for the displacements as functions of x, y, z 
and t. Since these functions are polynomial, and their derivatives 
are also polynomial, [f"] may be used in its "exact" form. 

Generally, the velocity-displacement law (2.31) will be linear. 


form of [g] is 

thus 

3/3. 

-03 


t 

Z 

y 

03 

3/3. 

-u 

z 

t 

X 

-03 

03 

3/a_ 


[g] = “z d/d t -“x (3.13 


where cu = (cj^, cu , 03 z ) is the angular velocity vector of the element 
as a rigid body with respect to an inertial reference frame, and, 
in general, is a function of time. Since this is linear, it is 
recommended that the inertial part of [K ] in (2.44a) be the same 
as that for a linearly elastic formulation. As a result, there would 
be no need for [N^] and {p}, and thus matrices [B^ ], [C^ G ], vectors 
(F^ } and (q^ }, and equations (2.45c,e) are eliminated from the analysis, 
Of course, in the future, one may wish to incorporate a true non-linear 
law (2.31). 


3.4 Constitutive Laws 

The bulk of the non-linear behavior will enter the formulation 
via the constitutive laws. The element under study will have the 
capability of analyzing a structure when it exhibits plastic or 
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creep behavior. In addition, parameters describing the constitutive 
laws will be temperature dependent. Linear behavior, in the form of 
elastic deformation and themoelastic expansion, will also be present. 


3.A.1 General Approach 

The general approach to the constitutive laws will be the 

assumption that the four forms of incremental deformation, i.e., elastic 

thermoelastic expansion, plastic and creep, are separable; thus, 

o o_ . , o_ , c> 

— = + -*P + -* C 

e e e e e 

Also, in many instances, it will be easier to express the constitutive 
laws using stress and strain deviator tensors, with pressure and 
dilatation measures as well; thus, the strain measure e, consisting 
of 7 quantites, is defined from e, the standard strain measures 
(3.10), by the relation 
e = [d]e 


where [d] is defined as 


[d] 


2/3 

-1/3 

-1/3 

0 

0 

0 

-1/3 

2/3 

-1/3 

0 

0 

0 

-1/3 

-1/3 

2/3 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

1 

1 

1 

0 

0 

0 


The new stress measure, S, is related to the old by the relation 


O = [d] T S 


(3.1^ 


(3.1 


(3.16 


(3.17; 


l n 
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It is thus seen that the energy integral term involving the inner 
product of stress and variation of strain is conserved, as is 
demonstrated by 

fitV - 6e T ([d] T S) - (6e T £d] T ) S * 6([d]e) T T-de T s 

o £ 

The constitutive law between o and e is written as in (2.34a). 
o o 

The law between S and e can be written as 

£ £ £ 
s = [E 7 ] e + 

where the 7 subscript indicates the size of the vector or matrix. 
Substituting (3.15) into (3.19) and noting (3.17) it is seen that 

X 6 = [d] T t 7 ; [Eg] = [d] T [E 7 ][d] 

where the 6 subscript, when removed, represents the conventional 
£ 

T vector and [E] matrix but is used here for emphasis. It will be 

£ 

seen that [E^] and are more easily manipulated than their 6-dimen- 
sional counterparts, thus, motivating the transformation. 


3.4.2 Elasticity 

It will generally be assumed that linear elasticity effects are 
always present, regardless of the other material effects that may 
be involved in the analysis. Linear elasticity is represented by 
specifying 2 of 3 parameters, E, G and/or V, which are related by 


G = 


E 

2(1+ v) 


The value of [E^] for linear elasticity is a diagonal matrix given as 


(3.1 


(3.1! 


(3.2C 


(3.21) 



42 


[E 7 ] = G 



2 


2 


2(H-\» 

3(1-2 V) 


(3.2 


Of course, if linear elasticity is the only effect present, the 

o o 

formulation of section § 2.3 applies. The quantity T (or Ty) is 
zero. 


3.4.3 Plasticity 

A model for plastic behavior was chosen that differs slightly 
from classical theoretical models of plasticity. First, uniaxial 
behavior of a virgin material is governed by a Ramberg-Osgood 
law (Ref [3.1]; specifically). 

The quantity Oy is the nominal yield stress of the material. Note 
that non-linear behavior begins immediately. The linear portion of 
the law reflects the elastic deformation. When the structure is 
unloaded, only the elastic deformation is recovered. Upon reloading, 
additional plastic deformation is assumed to occur immediately, the 
rate of which is determined by the stress state. See Figure 3.3. 


i +ji— 


(C 


(3.23) 
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Extension of this assumed uniaxial law to three dimensions is 
made by first recalling the definitions of e and S. Expressing 
these as tensor quantities, (3.15) and (3.17) become 

®ij “ £ ij “ 3 £ kk 6 ij 


c & kk 

S ij = °ij ” 1 °kk 6 ij 
p “ 3 °kk 

where e and p are the seventh elements of e and S respectively. 

The "effective stress" is defined as 

2 - 1/2 


a e “ t aS ij S ij + g P ^ 


A constraint is imposed on a g by stating that a g =a for uniaxial 
stress. Inserting a uniaxial state of stress in (3.25) yields 
6 = 9(1 - | a) 

3 

Note that gives the effective stress of 3^ theory. If the 

uniaxial rate law is given by 

oP c ,. .o 
e = f (o)o 


where for (3.23), 


f'(o) = j f (-2-) 
7 E 0 y 


n-1 


then its three dimensional extension is made by asserting 


3S 


ij 


I = f'(a ) a ^°e 
6 3p 


(3.2< 

(3.2^ 

(3.24 

(3.24 


(3.25 


(3.26) 


(3.27) 

(3.28) 


(3.29a] 

(3.29b) 


Noting (3.25) and (3.26) yields 
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9a e 

dS. . 

iJ 


as 


ii 


*°e _ aS kl S kl + 9(1 ~f a) P P 


3o 9(l-fa)p 

e * 5 

3p °e 

o 


3t 


<3. 3C 


Use of (3.30) in (3.29) gives 


M ‘ S Id ® kl + 

6 " {Sod-fa) P S kl s kl + 81(l-j a ) 2 p 2 ?) 


(3.31 


Generally, for an elastic-plastic element, (3.31) and the elastic 
constitutive relations are combined and inverted. In tensorial 
form, these relations are expressed as 

S ±j = 2Gej - y (ya 2 S ij + 9a(l-2/3a)S lj pi) 

p = 2G e - y (9a(l-2/3a)p e fcl + 81(l-2/3a) 2 p 2 e) 

y v 

where 

y = 3(l-2v ) 
l+o 

Y = f'(Op) 
e 

4 ' + 25 “ 2 Sfcl S,d + 15 ^ 81(1-2/30 2 p 2 

N y e e 

Note that the first terms in (3.32) are simply those that appear in the 

elasticity relations; thus, [Ey] may be expressed as 

[Ey] - [Ey] - B [eP] 

where [Ey] is given in (3.22), and [Ey]canbe shown to be-written as 

[fiP] = Ytv)(v} T 


(3.32) 


(3.33a 

(3.33b 

(3.33c 


(3.34) 


(3.35) 



46 


where 

^ a Sj 

4T a S 2 
* / y’ a 
a s 4 
'/v a S 5 
vF a S 6 
9(l-2/3a)S 7 /^T 

and 6 is a parameter which equals 1 or 0 depending on whether there is 

loading (§ e >0) or unloading (o e <0) at the point in question. 

The non-linear algorithm must be employed in order to solve 

structural problems incprporating plasticity. Generally, if only 

o_ 

elasticity and plasticity are taken into account, T = 0. Since it 
will also be assumed that structural inertia, if assumed significant, 
can be modeled linearly, the non-trivial matrices to be derived are 
[K e ], [A^], [B^], and [C^], and, if needed, [Xg] and [A^]. The 
additional set of shape functions to be incorporated into the analysis 
to evaluate some of these matrices is [N^]. These shape functions 
will be chosen to be the 28 functions based on the polynomials (3.5), 
i.e., the same spatial shape functions used for the displacements. 

These 28 incorporate those that would be present in an elastic analysis. 
The variable set {oq}, then, represent nodal values of the stress 
at (local) time zero. Since [N^] is not likely to change for prismatic 
discretizations, or from iteration to iteration, the [C^] matrices 
may be calculated immediately. In addition, [N^] and [N] may be used 
in (2.44b) (and (2.47b)) to calculate [A*] (and [A^]). The difficult 
aspect of the plastic element is to calculate [S]. 
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Suppose an iteration of the non-linear algorithm has been 
completed; thus, there are at hand "current" values of {u} and 
(o 0 ). From the latter, values of c g at spatial node i, designated 
are found at t=0 and t=T. In addition, from (3.30) and (3.32) 

§ e *, may be found at these times. There is now enough information 
to generate a cubic Hermitian polynomial in time to approximate the 
behavior ofc^ 1 . Taking a derivative gives a quadratic function of 
time describing . The roots of the equation § e **0, if they exist 
in the interval (0,T), indicate, for the upcoming iteration, an elastic- 
plastic boundary at node i. (There is, of course, a maximum of two 

such boundaries for any element.) In regions where 8 g i >0, plastic 

loading occurs and 6=1; where aj"<p, elastic unloading occurs and 6=0, 

e i 

or, [E ? ] = [E ? ]. Suppose a plastic interval occurs between t=t^ and 

i i i i 

where OCt^t^T; then [E^] t the i superscript referring to the 
i^ spatial node point, can be evaluated at t=t* and t^t*. A linear 
variation is now used in the interval so that 


^ ■ “to'Hp*) + c ** 3 (?r) 


In the elastic regions, of course, [E^]is[E^], a constant. Thus, 
in general, [E*] is a discontinuous function of time for an elastic- 
plastic element, though, through (3.37), an interpolated expression for 
it is now available. The [E*] matrix function of time is used to derive 
a matrix [S*] by 


[S 1 ] = / [E^][r] [N)dt' 

0 ' 

Note that [S*] consists of spatial dependence through [N], and 
temporal dependence through [E^] and [N]. Though [E*] is discontinuous, 
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[S 1 ] is not, though the functional form is discontinuous. Equation 
(3.38) may be evaluated at t«T to yield the matrix [S*]. 

The matrix [E*] represents the "contribution" from node i to the 

entire [E^3 matrix. The spatial dependence of [E^] is incorporated 

through the expression 
28 

[E 7 ] - l N* [E 7 ] 

' i=l 1 

where the N^'s are the 28 spatial shape functions produced from 
combinations of the polynomials of (3.5). If (3.39) is substituted 
into (2.38a), then the summation and N* terms may be taken outside 
the time integration to yield 
28 

[S] = Z N* [s 1 ] 
i =1 X 

If (3.40) is used in (2.44a) (and (2.47a)) it is seen that the time 

integration may be performed first, thus: 

28 T 

[K ] = Z f N* {/[N] T [r] T [S i ]dt} dV 
e i=l V x 0 

where the inertia terms have been ommitted. (The [X°] matrix is 
derived in a similar manner.) In (2.44d), the [B q ] matrix is calculated 
as 

28 

[B§] = I ( /[ N ]VdV)[S*] 

i=l V 0 x 1 


3*4.4 Creep 

The uniaxial law that will be used for creep behavior is 
£ C = KO m t n 

where K, m and n are material parameters. The rate constitutive 
law is found by taking a time derivative of (3.43) with the 


(3.39! 


(3.40) 


(3.41) 


(3.42) 


(3.43) 
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assumption (to be dropped later) that the stress field is kept 
constant; thus, 

gc - nKoV- 1 

Solving (3.43) for t and substituting in (3.43) yields 

JL ® 

o c nK 1 ^ 

? * 1 1 
n - 1 


Equation (3.45) is now assumed to hold in general. This is called a 
strain-hardening model of creep behavior, with e c being the accum- 
mulated creep strain. 

In order to extend the model to three dimensions, the parameter 

as in (3.25) is used with a=3/2. Another parameter, e c , representing 

a measure of accummulated creep strain, is defined as 

1 

r2 C C 1 2 

e c = *3 e kl e kl^ 

If only elasticity and creep are present to a significant degree in a 
problem, then 

e kl = ®kl " e kl = e kl " 2G S kl 


may be used in (3.46). If e c ) is defined as 

A £-1 

3 T .n n 
4> = 2 ^ g e 

i-1 

n 


then the constitutive law is 

e° c = d>s 
kl 9a kl 


(3.44 


(3.45 


(3.46) 


(3.47) 


(3.48) 


(3.49) 


Note that no rate quantity appears on the right-hand side of (3.49). 



50 


Adding the elastic rate law and inverting yields the constitutive law 

,e- 


f - [E®]? -4> [E*][I 6 ]S 


where 


zeroes 


(3.50 




(3.5i; 


zeroes 


Comparing (3.50) with (2.34a) indicates that 

£ „ 

T 7 - [E ? ][I 6 ]S (3.52) 

When employed in the non-linear algorithm, the development above 
implies the necessity, assuming that inertial effects, if significant, 
can be modeled linearly, of calculating [K g ], [A^J, [B^], [C ] and, 
possibly, [X e ] and [A^], {F®} and {q^}. First, it should be noticed that 

[S(t)] = [E®][f'][N(t)] - [E®][f'][N(0)] (3.53) 

because of the time-independence of [E_]; thus, [K ] is the elastic [K] 

/ 6 

described in § 2.3, (if [f'] is linear) minus a correction term. (The 

integral for [X^] is similar.) The matrices [A^], [A^] and [C^] are 

developed as in the elastic-plastic analysis of the last session, except 

with [S] simplified as in (3.53). The difficult part of the creep analysis, 

then, is to determine T and to be used in the {F } and (q } calculations. 

I X T 
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Expression (3.52) is analyzed comparable to the elastic-plastic modulus 
in the elastic-plastic analysis; that is, x will be approximated by an 


interpolation polynomial over the space time domain. For each of the 


§i, 


28 spatial node points, the function T (t) is introduced such that 


2 28 i £i 

T “ E t 

i=l A 


(3.54) 


Now S 3- may be interpolated using an existing expression derived in a previous 
iteration. The function <f>^ may be linearly interpolated by evaluating it 
at t=0 and t=T. While the evaluation of is straight-forward, the values 
of e c may be difficult to determine. For instance, if a linear elastic 
analysis is performed as an initial guess, e =0, particularly if (3.47) is 
used in (3.46); then, e £ should be calculated as 


e = K {fa m/n (Odt'} n 
C 0 6 


(3.55) 


where t here refers to global time. In particular, for a time t in the 


interval [t^, t^ + ^], a local relationship can be found by 

e = K { 5^ m/n dt' +./ o Wn dt'} 
c 0 e tj e 


/J \n 1 

■ Kf (f) + I 


s -'" • / o e m/n dt:'} n 


(3.56) 


When t=tj +1> the recurrence relation 


j+1 

e_ = K { 




(t) v 


j+1 / 

S o” /n dt') 
e 


n 


(3.57) 


is derived. If the integrand in the second term is modeled by '•assuming 


a linear variation of a , then (3.58) becomes 

e 



ol H 


52 


i 


e j+1 = K 
c 


!(?) 


( a e T -o e 0) ( 1 -hn /n) 


(a e T )' 


l+-lf 

n -(a e °) n JJ 


where T=t. .-t. and a and a are the values of a at the beginning and 

JTi J C 6 C 

end of the interval. 

The integrals for {F®} and {q^} are performed by first integrating 

over time with the time integral of *P, denoted as t 1 , or evaluating at 

t=T, and then taking the sum of general spatial integrals, one involving 

each . 

A 

3.4.5 Thermal Effects 

The effect of temperature on a structure is of importance in aircraft 
engines. Generally, the effect is two-fold; first, the material's desire 
to expand when heated may give rise to thermal stresses if the body is 
restrained against this expansion, or if portions of the body wish to 
expand more than others; second, the change in material parameters governing 
the various aspects of the constitutive relationships. 

The temperature field 0 is an assumed function of space and time that 
is given for any problem. It will also be assumed that this function is 
in polynomial form; "standard" temperature is defined as 0=0. 

The thermal expansion enters the general non-linear algorithm as a 
term; specifically. 


-Ea0/(l-2v) ■> 
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where the parameter a, the coefficient of expansion, is introduced. The 
dot over 6 represents the time derivative. If E, a and v are independent 
of time, then the time integration of (3.59) is straight-forward, and (3.59) 
gives rise to a {F®} and {q^} term. 

If the material parameters are assumed temperature dependent, then 
adjustments must be made accordingly. Specifically, the following parameters 
are chosen as functions of temperature: 


E - E 0 - E 1 e 

(Young's modulus) 

(3.60< 

“ * a 0 " “l 9 

(Coefficient of expansion) 

(3.601 

0 1 q 2 
’ °y - °y 9 

(Nominal yield stress) 

(3.60c 

K - K 0 /l 0 

(Creep parameter) 

(3.60d 

m - ra 0 e 

(Creep parameter) 

(3.60e 


if e<e 

a 

if 8_<0<0, (Creep parameter) (3.60f 

a o 

if 0 b <0 

where E Q , E^ a Q , c^, 0®, 0*, K Q , K^, m Q , ny n Q , and n 2 are 
material parameters. Parameters (3.60c-f) are used only in procedures 
that involve interpolation based on nodal point data, and thus, for 
(3.60d-f) in particular, their non-polynomial form is of Jittle 
concern. Equations (3.60a,b) may be substituted directly into expres- 
sions where E and a already exist in linear form. 


( ° 


n = 


n r n 2 
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TEXTRON 


3*5 Applied Load Vectors 

This section describes the various forms of mechanical loading 
that the element may undergo, and the analytical methods that may 
be used on the element level to find the proper contribution to 
{F*}. It must be remembered that {F*} is composed of a sum of terms, 
some of which are known and others unknown at the elemental level. 

This section addresses the known forces. 

3.5.1 Body Forces 

Body forces are forces per unit volume exerted over the entire 
element. They are represented by f. Generally, the body forces of 
interest in the applications that the 336-d.o.f. slave finite element 
may be used for are the weight per unit volume of the element, which 
is generally independent of time, and the inertial loading, i.e., 
the correction necessary to account for using a reference frame where 
the undeformed body is at. rest. 

The form of f may be either a specific polynomial in x, y, z 
and t, or tabular data at the node points from which an interpolation 
polynomial may be generated. The equivalent load vector, then, is 
generated from (2.13), where [N] is the shape function matrix for u. 

3.5.2 Surface Loads 

By surface loads, it is meant those force per unit area loadings 
that occur on the upper or lower surfaces only. Typical for these 
loads are fluid pressure or the placement of another object on the 
element. These loads will be represented by q, a vector function of 
x, y and t. Denoting the surface area by A, the applied load vector is 
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{F_} - / / [N] T q dA dt (3.61) 

b 0 A 

Like f, q may be a specific polynomial in x, y and t, or an inter- 
polated polynomial in these variables based on nodal data. The 
shape functions N are first evaluated at z** + 1/2 h(x,y), the sign, 
depending on whether the upper or lover surface is being used. 

3.5.3 Edge Loads 

By edge loads, it is meant those force per unit area loadings 

that occur on the surfaces that the elements may be connected to others, 

i.e., the 4 surfaces perpendicular to the x-y plane. These loads are 

usually due to the action of neighboring structures, and will be 

denoted by T, a vector function of £, the edge coordinate, z and t. 

Denoting the edge area of S, the applied load vector is 

T 

(F } = / / [N] T T dSdt (3.62) 

E OS 

The vector function T may be a specific or interpolated polynomial. 

When done this way, the [X^] and [A^] matrices should not include 

this edge as part of S^. An alternate method, which would leave the 

calculation of those matrices intact, would be to use the expression 

T 

{F } = / / [N r ] T T dSdt 

E 0 s B 


(3.63) 
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3.5.4 Point Loads 

If a load distribution is concentrated over a small volume or 

surface area, it generally can be modeled as a point load. If a 

particular point load acts at a (spatial) point which is not shared 

with any other spatial prism, then an equivalent load vector may be 

found at the element level. If the spatial point is (Xq, y^, z^) and 

the vector representation of the load is P(t) , then 

T 

(F p ) = / [N(x Q ,y 0 ,z 0 ,t)] T P“(t) dt 


(3.64) 


Generally, (x^, y^, z^) will be one of the 28 nodal locations, and 
thus [N(x 0 ,y^,ZQ, t) ] will have only 12 non-zero entries. 

If ( x o»yQ ,Z 0^ is s ^ ared with other elements, then (3.64) may 
only be used at the global level. 

The function P(t) may be a given or interpolated polynomial in 

time. 


3.5.5 Impulsive Load Distributions 

Impulsive load distributions may theoretically be caused by 
shock waves or a sudden jerking motion of the body. Generally, the 
first kind may manifest itself as a travelling pressure or point 
load. The second kind will generally be a body force load at a 
specific time. 

Figure 3.4(a) graphically displays the travelling pulse load. 
Since velocity is not continuous across the line in the space-time 
domain, the current element, using a prismatic discretization, could 
not be used in the neighborhood of the pulse. Figure 3.4(b) repre- 


May 1986 
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sents the second type of Impulsive loads. Though velocities are not 
continuous across this boundary, prismatic elements may be used. 

Hermitian elements may also be used if the precautions described at 
the end of § 2.5.2 are taken. The impulse per unit volume, J, has 
its load equivalent as 

£F t ) - f [Na] T J dV (3.65) 

J V 

where [N^] are the values of the shape functions over the volume 
at the time of the impulse. Equations (3.65) are only used on the 
global level. 

3.5.6 Point Impulses 

A hammer blow may be modeled as a point Impulse. Generally, 
a point impulse is the equivalent to a point force in conventional 
finite elements, and occupies the global right-hand side in a direct 
manner. As in conventional finite element analysis, it would be 
unusual if the point impulse was not applied at a (space-time) 
nodal point. Note that velocities would not be continuous at that 
nodal point. 

3.6 Integration Procedures 

Throughout this report, reference has been made to integrals over 
time and space. The spatial integrals can be over an edge, surface, 
or volume. In this section the procedures for taking these integrals 


are described. 
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3.6.1 Edge Integration 
First, it should be- noted that 

dS - dz dC (3.66) 

The z integration is performed first. The thickness along the 
edge is assumed linear; specifically, 

h(0 - h i (l- |) + hj | (3.67) 

where Si is the length from node i to node j . The limits on the z 
integfation are - 1/2 h(£) to + 1/2 h(£). The result is a new 
polynomial function of £, which is then integrated from 0 to 2. 


3.6.2 Area Integration 

The area integration is complicated because the area is assumed 
a general quadrilateral in the x-y plane. Consider Figure 3.5, where 
a quadrilateral has been divided into 2 triangles, I and II. For 
each triangle, a new coordinate system is formed with origins at the 
remaining corners. The £ coordinate (not related to an edge coordinate) 
is perpendicular to line 2-4, while the r| coordinate is parallel. Note 
that the directions of the and coordinates are opposite to 
their respective counterparts in triangle I. The integrands are 
polynomial functions of x and y (and, possibly t) . The linear trans- 
formations implied by the drawing indicate that for triangle I, the 
substitutions 

x “ + U I (y^-y 2 ) “ ti i (x 2 -x 4 )]/d 

y * y x + [C i (x 2 ~ x 4 ) + T i I (y 4 “y 2 ^/D 

be made, where D is the length of line 2-4. Defining 


(3.68) 
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n i 2 " C(x 2 - x i)(x a -x 2 ) + (y 2 -y 1 ) (y 4 -y 2 ) 3 /d 
\ ■ t( x 4“ x 1 )(x 4 - x 2 ) + (y4-yiHy 4 -y 2 )i/ D 

4 

5 I " r x 2 y 4 ’ x 4 y 2 + x l y 2 ” x 2 y l + x 4 y l “ X i y 4 ^ /D 

then the integral over the area of triangle 1 is performed by doing 

the n T integration first from n T S T /£ T to n_ £_/£•_ > resulting in 
l i 2 1 If 14 I If 

a polynomial function of and then integrating the resulting 


expression over from 0 to £ . 

For triangle II, use 

x * x 3 - C 5 II (y 4 -y 2 ) - t i ii (x 2 -x 4 )]/d 

y = y 3 " K II (x 2“ x 4 ) + n iI (y 4‘ y 2 )]/D 
n il2 = + (y 2 - y 3)( y 2- y 4^ /D 

n n 4 = t (x 4~ x 3 )(x 2'V + (y 4 -y 3 Hy 2 -y 4 )^ /D 

*II f - [x+y 2 - x 2 y 4 + x 2 y 3 - x 3 y 2 + x 3 y 4 - x^/D 

where the integration is first performed from to 

r) TT £ TT /£__ and the resulting polynomial in is integrated from 
112 J-J- llff *■ II 

0 t0 5 II f ‘ 


(3.6S 


(3.70 


3*6.3 Volume Integration 

The volume integration is essentially a z integration through 
the thickness followed by the area integration as described in the 
previous sub-section. In order to perform the z integration, the 
thickness h(x,y) must be specified. This is accomplished by inter- 
polating h separately over the two triangles I and II. This gives 
two separate polynomials h^. and h^, each linear in x and y, and 
also assures a linear h over each edge. The integration over the 
volume I is thus a z integration from - 1/2 h^(x,y) to + 1/2 hj(x,y) 
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and then the area integration over triangle I. The volume II 
integration is performed. in a similar manner, and then added to the 
volume I integration to yield the volume integral. 

As a closing note, it should be mentioned that the time Integra 
tion may be performed before or after the spatial integrations, 
except where specified. For non-prismatic elements, this would not 


be the case. 
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4.0 EXAMPLE PROBLEMS 

Sample problems were carried out for the example depicted In 
Figure 4.1 In order to demonstrate the capability of the element. First, 
a linear elastic transient analysis was carried out to validate the element. 

Two cases were examined. The first used a single element to cover the entire 
time, range of interest, while the second used two elements, each covering 
half the time interval. The results for the veritical displacement at the free 
corner are presented in Figure 4.2 and are compared with a NASTRAN transient 
analysis. The latter curve reveals that three inflection points are passed, 
thus making it difficult for even two cubic elements to model the behavior 
closely; yet, the results are more than reasonably accurate. 

Figure 4.3 demonstrates the convergence characteristics of the nonlinear 
algorithm for slightly plastic (o^ = 100 ksi) and highly plastic (o^ = 1 ksi) 
elastic-plastic constitutive laws for the vertical displacement at the free 
corder at the final time. Note that both curves begin with the linear 
elastic result for iteration 1 and oscillate about the final answer until 
convergence is met. These examples did not converge as fast as the truss 
examples presented in [4.1], though the slower convergence for the highly 
plastic case was expected. Figure 4.4 represents the final time history of 
the point loaded node for the two cases. Note that the oscillatory behavior 
demonstrated by the elastic case in Figure 4.2 does not appear to be present. 
This is probably due to the fact that the "instantaneous fundamental 
frequency" of the plate is so much lower dur to the reduction in the [E] 
quantities. 




p(t) = 10,000 t psi 

E = 3 x 10^ psi 

p = .1 lb/in^ 

v = .3 

o y = 1000/100,000 psi 

o =1. 

T = .001 sec 


o(0) = 0; u(0) = 0 


/ 


Figure 4.1. 


Sample Problem 
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5.0 CONCLUSIONS, SUMMARY, AND RECOMMENDATIONS 

In this section, an assessment of the successes and failures of the 
project will be made, as well as recommendations for future related work* 

The principal success of the project was the development of the temporal 
finite element method for nonlinear analysis. This new method presents an 
alternate approach to the standard "march in time" algorithms in use today. 
The major advantages are the ease in time step variation throughout the 
problem, progressively more accurate delineation of elastic vs. plastic 
loading in time, the natural way in which impulsive loading is handled, and 
the potential for true space-time coupling with the right choice of shape 
functions. While the convergence characteristics may be deemed disappointing 
it should be remembered that they must be compared with the rather small time 
steps used in most algorithms for nonlinear analysis. When more favorable 
meshing schemes are used, the convergence rate improves. 

There are many avenues available for future research using these methods. 
They include : 

• Improvements in the nonlinear algorithm. 

• Theoretical developments for characterizing initial value problems. 

• "Four-dimensional solid mechanics". 

• Implication and use of nonprismatic meshes. 

• Automated remeshing. 

• Applications to transient processes, including control. 

• Incorporation of the force method/hybrid method into the algorithm. 

• Expanded element library. 

A more moderate success was achieved in incorporation of the slave 
concept, i.e., the maximum use of the displacement shape functions and exact 
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Integrations. These conditions were followed as closely as possible in the 
derivation of all element matrix and vector quantities, and expressions were 
derived in closed form for each element of those quantities; however, the 
following combination of factors contributed to making the resulting 
expressions impractical, by nature of their length: 

• Four dimensions 

• Higher order shape functions 

• Nonlinearity 

• Nonrectangular geometry 

• Use of an older symbolic manipulator (FORMAC) 

When the geometry was restricted to rectangular the expressions were of a 
much more manageable length. The computer programmers working on the project 
believed that this was a more significant factor than replacing FORMAC with, 
say MACSYMA. 

Further research should be carried out in the slave method because of the 
tremendous potential accuracy improvements it may afford for a given set of 
shape functions. It is believed that the real fruitful work will be in 
improving and expanding the capabilities of the symbolic manipulator codes, 
rather than in the application. 

Whatever failures may have occurred in this program can be directly 
attributable to the size of the element and the scope of the characteristics 
it was to possess. It Is believed here that although this element is a 
suitable tool for the ultimate purpose of analyzing engine structures under 
severe thermo-mechanical environments, it was not a good candidate for the 
"proof -of -principle" work demanded by the developing technology produced under 
this contract. The recommendation here is to take a step backward to analyze 
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smaller, simpler elements in order to fully grasp the implications of the 
various methodologies. 
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